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ON  THE  ACCURACY  OF  A  NUMERICAL  INTEGRATION  PROCEDURE 
FOR  COMPUTING  DIFFRACTION  FIELDS  USING  TABLE-LOOK-UP 


INTRODUCTION 

In  a  classical  electromagnetic  problem,  it  is  often  required  to  compute  the  radiation  field  either 
radiated  directly  from  a  known  current  source  or  scattered  from  a  secondary  source.  In  the  latter  case, 
the  incident  field  on  the  scatterer  is  known  or  can  be  determined  with  the  known  boundary  conditions. 
The  vector  potential  function  of  such  an  electromagnetic  field  is  evaluated  by  integrating  the  current 
source  or  current  distribution  on  the  surface  of  the  scatterer.  Then  the  radiation  field  in  space  can  be 
determined  from  this  vector  potential  function.  In  general  this  leads  to  the  evaluation  of  the  integral 
[1,  2] 

E(0.  ♦)  -  Jt  J((.  i,)expU?(f.  V,  0))d(dii,  (1) 

where  J  is  the  current  density  excited  by  the  source  of  the  scatterer  surface,  and  £,t)  are  the  coordi¬ 
nates  of  the  current  source;  4>  and  0  are  the  coordinates  of  a  field  point.  It  can  be  shown  that  the 
above  vector  integration  can  be  converted  into  a  scalar  integration  function  such  that 

F(9.  *)  -  r>)cxp\jG((,  r,.  d>,  9)]dfdif.  (2) 

In  general,  the  current  distribution  and  the  geometry  of  the  source  or  scatterers  are  very  complicated. 
Except  for  a  very  few  simple  caes,  it  is  impossible  to  integrate  the  above  equation  into  a  closed  form. 
This  integration  is  usually  performed  numerically  on  a  digital  computer.  Although  it  is  straightforward 
to  program  such  an  integration,  the  required  computation  time  may  sometimes  be  extremely  lengthy. 
As  an  example,  if  we  assume  a  geometry  for  which  the  current  source  is  not  too  complicated,  a  modest 
100  mesh  points  for  each  of  the  f  and  tj  dimensions  may  be  required.  Then  it  requires  10*  points  of 
summation  for  this  integration  for  each  field  point  at  a  given  0  and  4>.  Now  suppose  that  a  100  by  100 
mesh  points  are  needed  to  map  the  entire  radiation  space.  Although  the  F((,  t>)  current  function  does 
not  need  to  be  repeated,  the  complex  phase  function  (the  Green's  function)  is  nevertheless  a  function 
of  both  source  coordinate  and  field  coordinate.  Thus,  computation  of  10*  points  is  required  in  this 
example.  At  each  point  a  computation  of  sine  and  cosine  functions  is  required.  In  general,  computa¬ 
tion  of  this  sine  and  cosine  is  a  time-consuming  process  for  a  digital  computer.  Most  fast  machines  can 
probably  compute  these  sine  and  cosine  functions  and  its  argument  within  perhaps  100  fis.  To  compute 
one  set  of  output  data  for  this  example,  a  central  processor  time  of  2.8  h  is  required  just  for  the  sines 
and  cosines.  Such  an  exercise  is  certainly  not  cheap. 

There  are  many  ways  to  combat  this  problem.  For  example,  Eq.  (2)  is  a  generalized  Fourier 
transformation.  Therefore,  the  fast  Fourier  transformation  (FFT)  may  be  applied.  Unfortunately,  the 
exponential  function  G(9,  ♦,  tj.  f)  in  general  has  no  linear  relation  between  source  coordinates  and 
field  point  coordinates.  In  general,  a  prerequisite  of  the  FFT  is  that  the  phase  function  G  computed  for 
any  field  point  (0,  4)  must  be  contained  in  a  set  that  is  finite  and  the  whole  set  is  generated  at  a  certain 
0,  ♦  point.  This  requirement  sometimes  puts  a  contraint  on  the  choice  of  the  integration  mesh  points 
that  cannot  be  achieved. 

There  are  other  ways.  One  example  is  found  in  the  computation  of  reflector  antenna  patterns.  In 
this  case  Gaiindo-Israel  and  Mittra  13]  proposed  to  expand  the  integral  with  an  infinite  series  as  a  func¬ 
tion  of  the  field  point  coordinates.  The  coefficients  of  such  a  series  are  the  Fourier  transformation  of 
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the  if,  i  coordinates.  To  compute  such  coefficients,  similar  integration  of  the  phase  function  is 
required.  Furthermore,  this  series  expension  method  applies  only  to  this  special  case. 

The  third  way,  perhaps  the  simplest  approach,  is  table-look-up.  In  this  approach,  at  the  beginning 
of  the  computer  program,  sine  and  cosine  functions  with  a  2v  period  are  computed  and  stored,  using 
an  adequate  number  of  sampling  points.  In  later  computations,  the  required  sine  and  cosine  functions 
are  found  from  the  table,  thus  avoiding  the  repeated  computation  of  such  functions.  However,  unless 
the  arguments  of  these  sine  and  cosine  functions  coincide  with  the  sampling  points  in  the  table,  errors 
will  be  introduced.  The  amount  of  error  is  a  function  of  the  number  of  sampling  points  in  this  sine  and 
cosine  table.  Therefore,  one  must  decide  how  much  error  one  can  tolerate  and  then  one  can  determine 
how  big  a  table  one  must  have.  In  this  report  we  analyze  this  problem  and  present  the  results. 

TABLE-LOOM-UP  ERROR 

For  a  numerical  integration,  Eq.  (2)  can  be  transformed  into  the  form 

Fi9.  <P)  -  p,  exp(/tjr„),  (3) 

It 

where  A„  -  f„)  >  0 

Gw  -  Gi-n,,,  9.  d). 

For  simplification,  the  double  summation  of  17  and  (  in  Eq.  (2)  is  converted  into  a  single  summa¬ 
tion  index.  The  function  F(9,  ♦)  is  a  complex  function  having  both  real  and  imaginary  components, 
thus 


F(9,  ♦)  -  x  +  Jy. 

When  table-look-up  is  used  to  determine  the  phase  of  the  exponential  term,  errors  are  introduced,  and 
Eq.  (2)  becomes 

F(9,  ♦)  -  £  A"  «xP(/G»)**P(/®»).  (4) 

tt 

where  the  8„  are  independent  random  variables.  Therefore,  An  exp(/G„)exp(/8a)  are  also  independent 
random  variables,  and  the  function  F(9,  ♦)  is  the  sum  of  many  independent  random  variables. 
According  to  the  central  limit  theorem  [4],  the  probability  density  function  of  Fi9,  ♦)  is  asymptotically 
normal.  Since  a  normal  distribution  is  completely  determined  by  the  first  and  second  moments,  our 
next  task  is  to  find  the  mean  and  the  variance  of  the  random  complex  function  F(9,  d). 

It  can  be  shown  that  the  mean  of  F(9,  d)  is 

fi9,*)-*il)'£iAHtxpijGn),  (5) 

It 

where  d(D  is  the  characteristic  function  of  the  random  variable  8W,  which  is  defined: 

d(*)  -  /  piS)exp(JkS)dS,  (6) 

and  pit)  is  the  probability  density  function  of  the  random  variable  8.  If  the  entries  of  the  sine  and 
cosine  table  are  sampled  uniformly  within  a  2  w  range  having  a  total  number  of  sample  points  A,  then 

piS)  —  1  fa  na  <  6  <  in  +  l)o  (7) 

—  0  otherwise, 

and  a  -  2  w/K. 
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There  are  two  ways  to  use  this  table.  The  first  is  to  choose  the  sine  and  cosine  function  of  the 
next  index  if  the  argument  of  the  phase  angle  is  larger  than  the  angle  na.  The  second  method  is  to 

choose  the  function  of  the  angle  na  if  the  argument  is  less  than  in  +  ~)a.  For  the  former  case  the 

characteristic  function  is 


*(*)  -  exp  -/y*a 


sin  y  fca 
2k* 


(8) 


while  for  the  latter  case  it  is 


*(*)- 


sinyfca 

lka 


(9) 


One  may  see  that  the  difference  between  Eqs.  (8)  and  (9)  involves  a  constant  phase  (1/2 )ka 
which  in  general  is  not  important  because  the  reference  phase  of  the  function  Fi0,Q)  is  arbitrary. 
However,  in  the  actual  programming,  the  first  method  is  somewhat  easier.  One  has 


2  *a  .  ,  1  (fa.)- 

~  "4  ‘ 


(10) 


Since  a  is  usually  very  small,  the  high-order  terms  are  ignored.  By  use  of  this  result,  Eq.  (5)  becomes 

F(9,  *)  -  (1  —  0.042a2)  £  A,  exp  «?,).  (11) 

n 

In  the  appendix  it  is  shown  that  the  variances  of  the  real  and  imaginary  components  of  FiO,  9)  are 
respectively: 

<r 2  -  }(1  -  *2(l))j;  A 2  +  -  *2(1)]£  A2  cos(2G„), 

and 

<r 2  -  Ul  -  *2(D) 1 4?  -  7^(2)  -  ^2(I)II  Ai  008(20,), 
while  the  covariance  of  xand  y  is 

o>  -  U+(2)  -  *2(1)J£  Ai  sin(2G„). 


The  joint  probability  density  function  of  the  real  and  imaginary  components  xand  y  is  then: 


pixy )- 


1 


2n<rx<r/J\  —  r2 


exp 


1 

(x-x)2  2 r(x-x)iy-y)  ,  (y-y)2 

2(1  —  r2) 

C r;  <rxay  <r; 

where  xand  y  are  respectively  the  real  and  imaginary  components  of  F{9,  4>),  and 


<r„  ■  r<r,<r 


xvy 


(12) 


(13) 


Our  goal  is  to  find  the  required  number  of  samples,  or  how  big  a  sine  and  cosine  table  is  required  to 
achieve  a  certain  prescribed  accuracy  for  both  real  and  imaginary  components,  x  and  y.  The  probability 
density  function  of  Eq.  (12)  is  too  complicated  to  estimate  this  relation.  However,  certain  approxima¬ 
tions  can  be  made. 


3 


I.  K.  HSIAO 


Since  £  A}  sin2<?,  and  £  >4*  cos2(/„  are  in  general  much  smaller  than  £ii,2  for  a  first  order  estima- 

mm  m 

tion,  one  may  assume  that 

<r2  -  <r  2  -  <r  *  -  y(l  -  *2(1))  £  >4  1 2.  (14) 


«rv*0. 

The  probability  density  then  becomes 


P (x  y)  -  exp 


I-p 


[(x  -  x)2  +  (y 


-  J)2)}- 


This  equation  implies  that  the  real  and  imaginary  components  of  F(9,  4>)  are  two  independent 
variables  and  have  the  same  variance.  Hence,  from  Eq.  (11) 

x  -  (1  -  0.042a2) Re  £  A„  exp (GJ,  (16) 

•nd  y  -  (1  -  0.042a 2)Im  £  A„  exp ((?„).  (17) 

m 

If  the  table  has  enough  sampling  points  and  a  is  small,  then  the  0.042a2  term  can  be  ignored. 

Both  x  and  y  are  then  equal  to  the  actual  value.  The  probability  density  function  for  either  x  or  y  can 

be  written  separately  as 


f»(x)- 


V2w<r 


(x-x)2 
2<r2  ’ 


when  a  is  small,  the  variance  can  be  approximated  by 


Let  us  assume  that  A,  is  normalized  such  that 


This  normalization  implies  that 


According  to  Cauchy's  inequality. 


£4,-1. 


I F(9,  *)|<1. 


Let  y„  -  1,  then 


H*  *  \MM 


Therefore, 
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Since  lx  -  x|  <  3 <r  with  a  probability  higher  than  99%,  if  one  desires  that  x  or  y  values  should  not 
deviate  from  the  actual  value  by  an  amount  c  with  a  probability  of  at  least  99%,  then  the  number  of 
sampling  points  of  the  sine  and  cosine  table  should  be 

*  >  7 W-  (23) 

where  A'' is  the  number  of  terms  summed  in  Eq.  (3).  For  the  example  given  in  the  introduction,  where 
N  “  100  x  100,  if  one  desires  that  the  error  e  <  10~6,  then  we  find  K  £  38  500.  Rather  than  com¬ 
pute  the  sine  and  cosine  terms  10*  times,  a  table  of  the  size  40  000  will  be  adequate  to  yield  results 
accurate  to  -120  dB  relative  to  the  absolute  sum  of  the  integrated  terms. 

NUMERICAL  EXAMPLE 

A  numerical  example  is  presented  here.  In  this  example,  the  antenna  pattern  of  a  uniformly 
excited  line  source  is  plotted.  Strictly  speaking,  this  is  not  quite  the  same  as  that  predicted  by  Eq.  (23) 
because  in  the  previous  analysis,  in  the  interest  of  general  scattering  problems,  the  scattered  field  is 
presented  as  a  complex  number  with  real  and  imaginary  components  separated.  In  an  antenna  pattern, 
the  amplitude  plotted  is  the  square  root  of  the  sum  of  the  squares  of  the  real  and  imaginary  com¬ 
ponents.  Probability  density  of  such  amplitude  is  known  as  Rician  distribution.  However,  the  illumina¬ 
tion  function  of  a  line  source  is  symmetrical  and  the  pattern  function  contains  only  the  real  component. 
Therefore,  the  results  of  Eq.  (23)  can  apply  to  this  case.  A  30-wavelength  line  source  pattern  is  plot¬ 
ted.  Since  this  pattern  function  is  a  sinx/x  function,  this  exact  pattern  function  is  used  in  both  Fig.  1 
and  Fig.  2  as  reference.  Next,  this  same  pattern  is  computed  by  numerical  integration.  A  tout  of  100 
points  is  used  in  the  integration.  Figure  1  shows  the  pattern  plotted  with  a  table  of  32  entries.  Accord¬ 
ing  to  Eq.  (23),  the  error  should  be  in  the  order  of  0.012.  In  Fig.  1  the  maximum  error  occurred  at  an 
angle  of  about  80°.  The  correct  pattern  has  a  —  36-dB  sidelobe  level  while  the  approximate  pattern  has 
a  peak  of  -31  dB  at  the  same  angle  point.  At  -36  dB  the  radiation  amplitude  is  0.0158  while  at  -31 
dB  it  is  0.0282.  The  deviation  of  the  amplitude  is  about  0.0124  which  is  very  close  to  what  Eq.  (23) 
predicts.  Figure  2  shows  the  same  antenna  pattern,  except  that  the  number  of  table  entries  increases  to 
128.  The  maximum  error  in  this  figure  is  about  0.0038  while  Eq.  (23)  predicts  0.003. 


CONCLUSION 

In  this  report  we  have  analyzed  the  error  introduced  by  table-look-up  for  the  numerical  integra¬ 
tion  of  an  electromagnetic  diffraction  integral.  We  have  shown  that  in  order  to  keep  the  error  of  both 
real  and  imaginary  components  to  be  less  than  c,  the  required  table  should  have  at  least  K  entries  uni¬ 
formly  distributed  within  a  2w  range,  and  K  >  3.85/cV^  where  N  is  the  number  of  terms  used  for 
numerical  integration. 
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Appendix 

COMPUTATION  OF  VARIANCE 

Let  F(B,  *)  -  x  +  Jy 

£{|F-£|2}-«r2  +  <r2 
E[(F  -  F)*\  -  <r \  -  &y  +  2 j<rv 
E[\F-  F|2}  -  £{|/|2J  -  |£(f)|2 

E[\F\2)  -  2  A 2  +  22  4,4.  expl/CG,  -  GJ]  •  l*(l)l2 

*  mm 

M*m 

l£(F)|2  -  |*<l)l2£  £4.4.  expire,  -  Gm)\ 

m  m 

*?  +  **-  (1-  1*0)  I2)  £4? 

H 

E[(F  -  £)2}  -  £{(f)2}  -  {£(f)|2 
£{(«}*- 24?*(2)expl/2(?J 

n 

+  22  4.4.  exp U(G„  +  G„)42(l) 

m  m 
n*m 

{£(/-))2  -  02(1)£  £  4.4.  exp\j(G„  +  Gm)\ 

m  m 

vx  -  ay  “  (4(2)  -  42(1)>2'<-  cos(2G„) 

It 

2<rv  -  (4(2)  -  420))24,  sin(2G„) 


